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Mischa Sallf] and Jan SmilQ 

Institute for Theoretical Physics, University of Amsterdam 
Valckenierstraat 65, 1018 XE Amsterdam, the Netherlands 
(Dated: February 1, 2008) 

The Hartree ensemble approximation is studied in the "symmetric phase" of 1 + 1 dimensional \(f> 4 
theory. In comparison with the "broken phase" studied previously, it is shown that the dynamical 
evolution of observables such as the particle distribution, energy exchange and auto-correlation 
functions, is substantially slower. Approximate thermalization is found only for relatively large 
energy densities and couplings. 
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I. INTRODUCTION 

Non-perturbative treatments of relativistic quantum 
fields out of equilibrium are currently under intense in- 
vestigation, because of their importance in applications 
to the physics of the early universe and to heavy-ion col- 
lisions. The Hartree approximation is sufficiently simple 
that it facilitates the study of inhomogeneous systems 
by numerical methods. The possibility of inhomogeneity 
allows for the incorporation of non-perturbative config- 
urations such as sphalerons and Skyrmions, or kinks in 
1+1 dimensions. 

Another aspect of inhomogeneity relates to equilibra- 
tion and thermalization. For homogeneous systems, the 
usual Hartree approximation fails to describe thermal- 
ization, because it does not incorporate direct scatter- 
ing. On the other hand, the classical approximation, in 
which expectation values are obtained as an average over 
realizations which are specified by an initial distribution, 
does allow for (classical) thermalization pj . The reason is 
that the realizations are typically inhomogeneous, even if 
the initial distribution is homogeneous, which allows for 
scattering of the classical waves. (Clearly, strictly homo- 
geneous classical scalar-field configurations cannot ther- 
malize as these correspond to a dynamical system with 
only one degree of freedom.) Classical waves that are 
small fluctuations on a ground state, correspond in the 
quantum analogue to a homogeneous mean field with a 
relatively small ( "vacuum size" ) two-point function. The 
classical scattering of such waves would correspond to 
direct scattering in the quantum analogue, which is ab- 
sent in the Hartree approximation. However, one may 
doubt, if strongly non-linear classical waves can be very 
well represented, in their quantum analogue, by two- 
point functions with, necessarily approximate, direct- 
scattering terms in their evolution equations. One can 
furthermore imagine that in such situations strong non- 



linearity is important for the first stages of equilibration, 
and that this aspect could be better represented by strong 
and typically inhomogeneous mean fields. 

Motivated by the capability of the classical description 
to incorporate strong non-linearity, we recently intro- 
duced ytyf a Hartree-ensemble approximation, in which 
the initial density matrix is expanded on a basis of Gaus- 
sian coherent states, 



D(j) Dtt p[(f>, 7r] |</>7t)(07t|, 



(1) 



with "mean fields" <f> and tt = <f> and suitable variances. 
With "mean field" we here denote the expectation value 
in the Gaussian quantum density matrix \<p, ir}(4>, tt\. The 
true mean fields are obtained after taking also the clas- 
sical average over the separate "realizations", weighted 
with the classical density functional p[<fi,Tr]. For each of 
these "realizations" separately we use the Hartree ap- 
proximation. 

This formulation achieves two things. First, it allows 
for the description of systems with non-Gaussian initial 
density matrices, thereby going beyond the Hartree ap- 
proximation, and second, it allows for strongly non-linear 
equilibration. There is also a remnant of the "small fluc- 
tuation type" scattering: the "mean fields" in the "real- 
izations" \4>tt) (4>tt\ are in general inhomogeneous, even if 
the full expectation values describe a homogeneous sys- 
tem, and these inhomogeneities lead to indirect scatter- 
ing of the quantum modes via the Fourier modes of the 
inhomogeneous "mean field" . 1 

Note that when the initial distribution p[</>,7r| is pos- 
itive, the initial correlation functions can in principle 
be calculated to arbitrary precision using Monte Carlo 
methods. However, the dynamical equations of motion 
still apply to an individual Gaussian "realization", us- 
ing a Gaussian approximation (Hartree), and the corre- 
lation functions at later times depend on these approx- 
imations. When the "mean field" dynamics in the "re- 
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Of course, we would like to include direct scattering as well, as 
in 0, but this is at present numerically too demanding in the 
inhomogeneous case. 
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alizations" dominates over the Hartree corrections, we 
expect that this Hartree-ensemble method works better 
than the usual Hartree method in terms of the true mean 
field, especially in case the latter is homogeneous. 

In Ref. 0, 0| we studied this extension of the usual 
Hartree approximation in the 1+1 dimensional </> 4 model, 
and found that it leads to approximate thcrmalization. 
We focused on the "broken phase" of the system in these 
papers. In the present work we discuss the "symmetric 
phase" ? 

The reason for presenting results also in the "symmet- 
ric phase" is twofold. First, it is interesting in its own 
right. Not much is known yet about inhomogeneous rel- 
ativistic systems out of equilibrium and numerical data 
can help to sharpen our intuition. Second, results by Bet- 
tencourt et al. |5| , which were obtained in the "symmetric 
phase" , appear to contradict our earlier results in 0, • 
We try to clarify the situation here, and substantiate 
our earlier observation that the evolution towards equi- 
librium in the "symmetric phase" is much slower than in 
the "broken phase" . 

Broadly speaking, our results are as follows. We find 
(relatively) rapid approximate thermalization from an 
initial homogeneous non-equilibrium density matrix, in 
case of high energy-densities. We interpret this as con- 
cording our intuition on strong "mean fields" . On the 
other hand, for low energy densities and/or couplings the 
time scales for reasonable changes in particle distribution 
functions are huge, even larger than what we found ear- 
lier in the "broken phase" , which makes numerical study 
very difficult. The inhomogeneous case studied by Bet- 
tencourt et al. (a single Gaussian wave packet) falls in 
the low energy density class and indeed, even after very 
large times, we do not find thermalization. A more com- 
plete summary of the many detailed results will be given 
at the end of this paper. 

In Sec. |H] we will briefly review the equations of mo- 
tion, explain the Hartree ensemble approximation, de- 
scribe various initial conditions and observables, and de- 
scribe some aspects of the elegant definition of the theory 
on a space-time lattice. In Sec. lIIII we will present results 
obtained using numerical simulations. We give an expla- 
nation of a phenomenon found in the simulations that 
we call "local fc-space equilibration" in Sec. lIVI where we 
also comment on the difference between the "symmetric" 
and "broken phase" . A discussion is in Sec. [V] Finally 
in the Appendix we derive the initial as well as free-field 
form of the particle distribution for the Gaussian wave 
packet. 

In the rest of this paper we will drop the quotes around 



"mean field" . 



II. HARTREE APPROXIMATIONS, INITIAL 
CONDITIONS AND OBSERVABLES 



A. Hartree approximation 

For the \<p theory in 1+1 dimensions the Heisenberg 
equations of motion are given by 



{-d 2 + ^{x) + Xcj){xf 



0. 



(2) 



In the Hartree approximation (H.a.) we can write 4>(x) 
and its time derivative tt(x) as 



4>{x) H =' <j){x 



)+Y.[ h ^f^) + b a r a {x) 

a 

n(x) H = a ' n(x) + £ [bU a (x] + bj* a (x) 



(3a) 
(3b) 



where tfi(x) = (<fi(x)} and f a (x) are the time-dependent 
mean field and mode functions, and b a are time-mdepen- 
dent creation and annihilation operators satisfying the 
usual commutation relations. Without an approxima- 
tion a similar expansion is possible but only with time- 
dependent creation and annihilation operators: in terms 
of the initial b a and w a the dependence is then non-linear 
at later times. The fact that the same set can be used at 
all times is due to the Gaussianity of the Hartree approx- 
imation (see e.g. Q for more information). Combining 
Eqs. and © and taking expectation values gives the 
equations of motion for <fi and the f a , 

4> = A(j) -[/i 2 + A0 2 +3AC]</>, (4a) 
f a = Af a - [m 2 + 3A0 2 + 3XC]f a , (4b) 



with 



C = V(2< + l)|/ Q | 2 , n° a = (b a b a ). (4c) 



The exact equation of motion for the one-point func- 
tion, which follows by taking the expectation value of the 
Heisenberg equation of motion (J2J, contains the three- 
point function. An exact equation for the three-point 
function can be obtained from taking the expectation 
value of the product of Eq. J2J) with two field operators 
and will contain the five-point function. The resulting in- 
finite hierarchy of higher point functions is truncated by 
the Hartree approximation which factorizes these higher- 
point functions into one- and two-point functions. 



2 At zero temperature the expectation value of the order field (f> 
distinguishes the phase regions in coupling-parameter space. At 
finite temperature this expectation value vanishes due to kink- 
antikink condensation. However, there is still a clear distinction 
in the behaviour of the system in the two phases even in finite 
volume. 



B. Choice of mode functions and renormalization 

As in our previous work, we will use the stationary 
solutions of Eqs. J3J to define the initial mode functions, 
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i.e. we will choose them as plane waves with wave vector 
k (a — > k) 



AM) 



fk(x,0) = -iuj y k 



(o) 



( 5 ) 

Here L is the "volume" of our one dimensional sys- 
tem (using periodic boundary conditions) and uo k ^ — 
\J m? + k 2 , with m the physical mass, defined self- 
consistently as the square root of the term in square 
brackets in Eq. J4bj. The sum l(lcjl for these initial 
mode functions diverges logarithmically in 1+1 D, which 
is compensated by the bare p 2 . In the vacuum, i.e. us- 
ing stationary solutions for 4> and f a of Eqs. (QJ with 
nP a = 0, the renormalized mass parameter and mode sum 
are given by 

1 3AC„o - 



2 



(i 



/i 2 + 3AC* 



2 

Mrcn 



3AC r< 



(6) 



Then m 2 = p 2 en + 3Aw 2 , with v the stationary value of <f>. 

In this paper we will use n° a — 0. Using the mode func- 
tions (|SJ), the coherent states \cj>, it) satisfying 7r) = 
0, have mean fields 4>(x) and 7r(x), with a variance that 
is easily computed. We shall use these states as basis for 
our Hartree ensembles 



C. Initial conditions 

We used two different initial conditions for the mean 
field, a sum of standing waves with a flat distribution of 
phases, which we also used in our previous work in the 
"broken phase" Q, and a single Gaussian wave packet, 
as studied by Bettencourt et al. pj. 

The first is given by 



'(x) = 0, ir(x) = 



Jmax 

Am cos(2irjx / L — 

.7 = 1 



V>i), (7) 



where the maximum momentum 27rj max /L is typically 
of the order of the mass m and the constants ipj are 
random phases with a flat distribution (i.e. they are uni- 
formly distributed in [0, 2ir)). We shall call such p[(p, it] 
flat ensembles. The energy, which is independent of the 
phases ipj, is given by 



E _ A 2 Lmj n 
m 4 



(8) 



We use both A and j max to vary the total energy density. 
The second initial condition is a Gaussian wave packet: 



n(x) = 0, 

Its energy is given by 
E $ 2 



(x) = $ exp 



x 

2A 



- = ^J-^(2 + AAm 2 + V2AX^ 
™ 8 V Am 2 K 



(9) 



(10) 



We will restrict ourselves to Am 2 = 2 and use $ to vary 
the total energy in the system. For this type of initial 
conditions we do not average over multiple runs, so p[4>, it] 
is a delta functional and p is a coherent pure state. 



D. Observables 

As in we define particle numbers and frequencies 
LOk in terms of equal time correlation functions, 



S(x,y) = (<t>(x)ct>(y)) - (<t>(x)) (<%)}, 



U(x,y) = (fc(x)iT(y)) ~ (fc(x)) (w(y)). 



(11a) 
(lib) 



The over-line indicates coarse graining over all space and, 
depending on the simulation, also includes course grain- 
ing over a time interval and an ensemble of initial condi- 
tions. The symmetrized 7r0-correlation function, T(x, y), 
vanishes in equilibrium, and we will just use S and U . 
Taking the average over all of space leaves S and U de- 
pending only on the coordinate difference x — y. In terms 
of the Fourier transform, 



S k 



dx e~ ikx S{x, 0) 



(12) 



and similar for U, the time-dependent particle numbers 
and frequencies are defined by 



S k 



1 



(13) 



Therefore, apart from numerical corrections discussed in 
Sec. Ill El we define instantaneous particle number and 
frequency by 



VUkSk 



(14) 



In practice n k is positive (it can be shown to be posi- 
tive provided the symmetric correlation between <p and 
7r vanishes) . Other definitions of particle numbers are in 
use, e.g. derived from time-dependent creation and an- 
nihilation operators defined in terms of adiabatic mode 
functions, 

a fc (t)e- i ^'' i »( f ') = 

i r L 

dxe~ lkx [Cj k {t)(j>{t,x) +iw(t,x)]. (15) 



y/2Q k (t)L Jo 

However, one then still has to choose This can 

introduce some ambiguity, especially if the system is far 
from equilibrium and the effective mass term in the equa- 
tions of motion for the modes becomes negative; see for 
example [EH- Using Eq. ijTil) as the definition of particle 
number, has the advantage that the frequencies are real 
and positive by construction. 

The correlation functions S and U may be written as 
a sum of contributions from the mean field and mode 
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functions separately. Given ojk as obtained from the total 
S and U, we can study their separate contributions to 
as defined by Eqs. fT3|) . 

The energy density can be obtained from the effective 
Hamiltonian H c h = (H) |2j. We split it into mean field 
and mode contributions according to 



E mi =i(5 t 0) 2 + \(d x 4>Y + \& n <t> 2 + 

1 2 2 \ 4 



(16a) 



(16b) 



+ ■^^Cron0 2 + T^Crcn ~~ -^vac 



where v, the ground-state value of <fi, is zero in the "sym- 
metric phase" . Although the splitting between modes 
and mean field is somewhat arbitrary, the above defini- 
tion has the advantage that it reduces to the classical 
expression when the modes are of the vacuum form JSJ . 

The quasi-particle energy is also an interesting observ- 
able. It may be defined as 



k 



(17) 



where n k can be obtained from the two-point functions 
of the mean field, of the mode functions, or of the sum 
of both. 



where p and q are canonical operators satisfying [q, p] = i. 
A finite time evolution then takes the "Trotter form" 

U q U r =U q --- U p U q U p U q U p U q ■■■U q , (22) 

The Heisenberg operators 

p(t) = U r 1pU r , q(t) = U r] qtj r ', t = a r, (23) 

satisfy the discretized equations of motion in leapfrog 
fashion, 

p(t + ao) = p(t) - a V (?(*)), (24a) 
q(t + a a ) = q{t)+a p(t + a ). (24b) 

With q(t) — > q(t), p(t) — ► [q(t) — q(t — a )]/ao, the above 
Eqs. 1241) are identical in form to the classical equations 
obtained from the stationary action principle. 
Making a unitary transformation 



f _ e -ia V(q)/2 jj e ta„V(q)/2 



(25) 



we get an equivalent operator T, which becomes the 
Hermitian and positive transfer operator upon analyti- 
cally continuing to imaginary time (see e.g. writing 
a = e- iB \ao\, 8 = 0^ir/2, 

f _> e -\ao\V(q)/2 e ~\a \p 2 /2 e -\a \V{q)/2 _ ^6) 

Specializing to the harmonic case V(q) = uj 2 q 2 /2 we 
can diagonalize the time evolution in terms of creation 
and annihilation operators & and c, 



E. Implementation on a lattice 

The discretization of the scalar field theory on a space- 
time lattice has some elegant features which we present 
briefly in this section; for fermions, see 8]. For simplicity 
we start with a simple quantum mechanical system of 
unit mass, with action 



S 



a )-q(t)f 



2a 2 , 



V{q{t)) 



(18) 



where ao is the time step, t = aor, with integer r. We 
define the quantum system by means of the path integral. 
The discretized path integral 



Z 



[Jjdg(*) 



JS 



(19) 



corresponds to an evolution operator in Hilbert space 
that is a product of single step evolution operators given 
by 



with 



U p = e 



u = u p u q , 



■ia V(q) 



(20) 



(21) 



f cf f = e iaoUJ c, fctft = e~' laoUJ c\ 



with 



\/2w(") 

and 

cos(a Lu ( - e ' ) ) = 1 — idQW 2 , 



(27) 



(28) 



(29a) 



= J- s in(a w (e) ) = wWl - i a 2 w 2 , (29b) 
ao V 4 

and the conjugate relation for c'. The creation and anni- 
hilation operators satisfy the standard commutation re- 
lation [c, c'] = 1. The superscripts e and n distinguish 
the "exponent omega" (eigenvalue omega) a/ 6 ' from the 
"normalization omega" (eigenvector omega) w^ n \ and 
both go over to the "original omega" uj in the contin- 
uous time limit ao — > 0. The ground state is given by 



- ,„-«- f »V/a 



(30) 



c\0) = 0, (g|0) = ue 
with v a normalization constant and 

f(c t )™|0) = e -*(»+i/2)ao^) (ct)»|0). (31) 
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The evolution becomes unstable when a^uu 2 > 4, for 
which tu k is imaginary. The eigenvalues of T are then 
no longer phase factors and its eigenfunctions no longer 
normalizable, despite its formally unitary form. This is 
of course avoided by taking ao sufficiently small. The 
discretization errors in a/ 6 ' and lj^" 1 ' are of order a 2 ,. 

It is natural to identify the Hamiltonian H from T = 
exp(— iaoH), but this leaves a modulo 27r/ao ambiguity 
for the eigenvalues of H (the imaginary time version is 
unambiguous). To pin down H more precisely we can use 
the Baker-Campbell-Hausdorff series for combining the 
exponents in T, which gives H — p 2 /2 + V(q) + 0(o§). 
We shall neglect the corrections of order a 2 ,. The exact H 
is time independent. In practice, the expectation value 
of the approximate H is constant in time up to small 
fluctuations, as expected for a leapfrog algorithm. 

For application to the Hartree approximation it will be 
more convenient for us to work with the unitarily related 
creation and annihilation operators that diagonalize U, 



e ia V(q)/2 £ e -ia V(q)/2 

1 (l- e-' ta 



ia 



-q + ip 



(32) 



for V = oj 2 q 2 /2. Note that a — > c in the limit ao — * 0. 

The generalization of the above quantum mechanical 
model to our scalar field is straightforward. The lattice 
action on a space-time lattice with spatial-temporal lat- 
tice distance a/ao is given by 



S\d>] = 



X,t 



[<j)(x,t + a Q ) - (f>(x,t)] 2 
2a 2 

[(/)(x + a,t) - 4>{ x,t)f 1 2 2 

U (p(X,t) 



2a 2 



-\(j){x,ty 



(34) 



where we assume a periodic physical size L = Na. The 
operator description in Hilbert space follows from the 
lattice regularized path integral. In the Hartree approxi- 
mation we write the operator fields in terms of a complete 
set of mode functions, 

4>{x,t) = <p(x,t) + J2$kfk(x,t) +blf k (x,t)*], (35a) 
ft 

n(x,t) =tt(x, t)+ J2[hfk(x,t) + blf k (x,t)*}, (35b) 

k 

where the use of 

f{x,t) = (36) 

a 

is inspired by Eq. I|24b|) . (Using instead the forward 
derivative f k (x, t) = [f(x,t+ao) — f(x,t)]/ao gives equiv- 
alent results.) Imposing canonical commutation relations 



for both (j>, 7r and b k , b k , leads to the orthonormality and 
completeness relations 

aY%f k (x,t)ft(x,t) -i/fc(i, *)/,*(!,*)] = (37a) 

X 

~ iMx,t)f* k (y,t)} = 5 -f. (37b) 

* — • n 



The time independence of the orthonormality conditions 
corresponds to Noether charges of symmetries of the 
effective action on the lattice. We use the static solutions 
of the Hartree equations in constructing the set of mode 
functions. Their equation of motion 

f k (x,t + op) - 2f k (x,t) + f k (x,t- op) _ 
a 2 

fk(x + a,t) - 2f k (x,t) + f k {x - a,t) 2 
= m f k [x,t), 

(38) 

can be written in the leapfrog form l|24|l . The solution of 
the recursion relation l|38(l can be written as 



(33) fk(x,t) 



„ikx — iuj) t 



2J k n) L 



2irj N N 

k ~—> 3 ~ ~T + 1 '" ' ' T' 

(39) 



giving 



2 — 2cos(w^ao) 2 — 2cos(fca) 2 



+ m 2 = 0. (40) 



Defining a lattice ui^ as 



(a) 



2 — 2 cos(fca) 



(41) 



we find the analogue of Eq. i|29a| 



cos(a ^ e) ) = l-i^(4 Q) ) 2 , (42) 

which has real (Jjf 1 solutions for a/ao > \/4 + a 2 m 2 . We 
used a/ao = 10 in our simulations, which amply secured 
stability. The normalization in Eq. Ij39(l is fixed by the 
orthonormality relation (|37a|l . which gives the analogue 
of Eq. pjbj 



(n) _ sin(a 4 e) ) _ (a) f 



(43) 



The completeness relation (|37bll is then also satisfied. 
When the mode functions have the form Ij39(l . the a k 
defined by $ = J^ k a kfk + h.c, tt = J2 k " a *.h + h.c, are 
related to </> and tc as in the quantum mechanical case 



Note that 



w fc ^k 



oj), in the limit ao — > 0, 



and wl a ^ — > V m 2 + fc 2 as a — » 0. 
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We end this section with a properly discretized version 
of the instantaneous particle number rife, using the sta- 
tionary solution and the two-point functions (fTTJ) . 
Suppose the mean field is zero and 



(Kbk) = nl = n°_ k . 



Then 



Sfc(t) 



U k (t) 



where we used 



f k (x,t)fUv,t) = (4 a) ) 2 f k (x,t)r k (y,t). (46) 

Inverting Eqs. 1)45(1 we find that our definition of the 
instantaneous energy uikit) does not need discretization 
corrections, 



in)' 



K 



(a)\2 
k ) 



,(«) 



(44) 

(45a) 
(45b) 



,(«) 



' u k (t) 
s k (t) 



ui k {t). 



(47) 



On the other hand, the definition of instantaneous par- 
ticle number needs important corrections for large u>f.' 



yUkSk 



(n) 



,(<*) 



Uk(S k ~ -alU k ) = n k {t) + 



1 

2' 
(48) 

using Eqs. lO and 63 • 

For larger energies the corrections can become quite 
important. Denoting the uncorrected particle number 
by n k = y/UkSk - 1/2, we find 

h k - n k n k + ' 



n k 



nfc + I 1 ..(aha 



(49) 



Using a Bose-Einstein distribution at T = to and the typ- 
ical value qqttl = 1/80 we find that the relative difference 
becomes unity for uik/m = 7.5. At the lower temperature 
T/m = 0.5 this is the case already for oj k /m — 4.3. 



III. NUMERICAL RESULTS 

In this section we will present data we obtained using 
numerical simulations in the "symmetric phase" , i.e. u 2 
and A corresponding to the "symmetric phase" at zero 
temperature. First we will discuss the particle distribu- 
tion to study its equilibration behaviour and to search 
for thermalization. Then we will examine the energies 
and auto-correlation function to analyse the time scales 
in the theory. 



N=128, Lm=32, Wm 2 =1/6, E/Lm 2 =4, Modes only, 20 initial cond. 



- P(o-M-) 




tm=0 




tm=1000 




- tm=10000 




pm=0.41 




" |i/m=0.59 







0.5 1 



1.5 



2.5 



3.5 4 0D k -> 



FIG. 1: The particle numbers log(l + l/n k ) in the modes 
as a function of u>h- Average over 20 flat ensemble initial 
conditions, A/m 2 = 1/6, E/Lm 2 = 4, at times up to tm = 
10 4 . Time increases from the top curve to the bottom curve. 



A. Flat ensemble 

In the flat ensemble of initial conditions, the initial 
mean field <f> of a realization is equal to its vacuum ex- 
pectation value 0, while its momentum is the sum of 
waves with random phase, as specified in Eqs. JJJ. We 
averaged over 10 or 20 initial conditions, and excited all 
non-zero modes up to /c max = 2irj max ,/L = Trm/2. Sim- 
ulations have been carried out for three different cou- 
plings A/m 2 = 1/6,1/8,1/12 and three different energy 
densities E/Lm 2 — 4, 2, 1, as well as for the combination 
A/m 2 = 0.1 and E/Lm 2 = 0.4. We mainly used N = 128 
lattice points, and volume Lm = 32, while the temporal 
lattice distance ao = a/10. 



1. Particle distribution function 

Fig-fflshows the particle number obtained for coupling 
A/m 2 = 1/6 and energy density E/Lm 2 — 4. As in our 
previous work we compare the out-of-equilibrium particle 
densities with a Bose-Einstein (BE) distribution 



1 



nk 



(50) 



We therefore plotted log(l + l/n k ) versus uj k /m, since 
a BE distribution then shows up as a straight line with 
slope m/T and offset —a/m (T temperature, \i chemical 
potential). For this largest coupling and energy density 
in our study we find approximate thermalization to the 
BE form with a temperature T/m — 2.4 and chemical 
potential u/m = 0.6. In contrast to what was found 
in the "broken phase" 0> Q , a substantial chemical po- 
tential is needed to make a reasonable fit. Another dif- 
ference is the larger time scale involved: in the "broken 
phase" at an energy density E/Lm? = 0.5 and the same 
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E/Lm 2 = 1 


E/Lm 2 = 2 


E/Lm 2 = 4 


A/m 2 


= 1/6 


j3m 


1.12 


0.71 


0.41 






fi/m 


0.95 


0.83 


0.59 


A/m 2 


= 1/8 


j3m 


0.89* 


0.60 


0.45 






fi/m 


0.62 


0.76 


0.80 


A/m 2 


= 1/12 


j3m 




0.68* 


0.40 






fi/m 




0.76 


0.85 



TABLE I: Inverse temperature j3 and chemical potential fi 
as derived from a Bose-Einstein fit to the particle numbers 
(modes only). See text for further explanation. 



VlMrenl = (A/m 2 = 1/12), we could already recog- 
nize BE behaviour with T/m w 1 at a time tm < 100, 
while here, at an 8 times larger energy and roughly 2 
times larger effective BE temperature we can only clearly 
do so at time tm > 2000. A fit of the local tempera- 
ture T(t) approaching approximate equilibrium gives an 
equilibration-time scale totbe = 1500—1600 (exponential 
fit over k/m < 1.7, 100 < t < 6000). 

For most parameters used in our simulations we can 
recognize Bose-Einstein features in the low momentum 
part of the distributions, and linear fits can be made as in 
Fig-H The results of these fits are shown in TableQ] Fits 
marked with a star were made at tm = 59000 • • • 60000, 
all others at tm — 9000 ■ • ■ 10000. Including the mean 
field in the two-point functions, typically gives the same 
temperature within errors, but a noticeably larger value 
(+ 0.05) for jit/m, corresponding to a higher particle num- 
ber. The chemical potential is also more sensitive to the 
exact fit-interval than the temperature. In the two *- 
marked runs a thermal distribution could be recognized 
only at times > 20000 (A/m 2 = 1/12) and > 45000 
(A/m 2 = 1/8). 

For the run at A/m 2 = 1/12, E/Lm 2 — 1, we did not 
find a thermal-like distribution even at the latest simu- 
lation time tm = 10 5 . We see that the energy is trans- 
ferred from the mean field to the modes and the system 
equilibrates "locally in fc", but the total particle num- 
ber remains roughly unchanged. The same was found 
at the lower energy density E/Lm 2 = 0.4, and also for 
the Gaussian wave packet initial condition. We interpret 
this as a resonance phenomenon in the equation of mo- 
tion of the mode functions, which will be described in 
Section ITVBl 

Comparing the results at E/Lm 2 = 2 and 4 it seems 
that the temperature only depends on the energy density 
and not on the coupling. This appears to hold even for 
the distribution function itself, as illustrated in Fig. |2 
where the distributions for different couplings are plot- 
ted at different times. The different times at which the 
curves in the figure overlap suggest that the equilibra- 
tion time scale for the particle distribution is propor- 
tional to A~ 3 . The same power is found at the energy 
density E/Lm 2 = 2. Table [I] shows that the temperature 
is roughly proportional to W E / L, which can be under- 
stood from the scaling behaviour in Fig. [21 since there is 
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FIG. 2: Scaling behaviour of the particle distribution at fixed 
energy density and differing couplings and times. 
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FIG. 3: The particle numbers log(l + l/nfc) in the modes as a 
function of ujk at large times tm > 10 4 . In the small k region 
time increases from bottom to top, in the large k region time 
increases from top to bottom. 

no other scale left. The same argument should apply to 
the chemical potential. However, this quantity is more 
dependent on time than the temperature and runs at dif- 
ferent parameters are best compared at different times as 
in Fig. El which we have not done in Tabled 

The independence of coupling suggests that a represen- 
tation of the energy in terms of almost free quasi-particles 
will be reasonably good. We will check this in the next 
section. 

Fig. [31 shows the distribution for late times, when it 
starts to deviate from the BE form. Note the difference in 
vertical scale compared to Fig.[l] At tm = 15000-20000 
classical-like deviations become visible in the form of con- 
cave behaviour at low u>k {nt = T/uk => (d/dujk) 2 log(l+ 
l/n/0 < 0). 

The mean field in this time region behaves very inter- 
estingly. In Fig. 01 we plotted the particle numbers at 
tm = 50000 ■ • • 70000 as a function of momentum fc, both 
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N=128, l_m=32, Wm 2 =1/6, E/Lm 2 =4, tm=5-10 4 ...7-10 4 , 10 init. cond. * N=128, Lm=32, X/m 2 =1/6, E/Lm 2 =4, 10 initial cond. 




FIG. 4: Log-log plot of the particle numbers for mean field 
alone, and combined with modes, vs k/m, at large times. 



FIG. 5: Contributions to the energy density, for short times. 
^,From top to bottom at tm = 450: JZfc nfea -' fc (mean field 
+ modes), E from H eS , J] fc n fc w fc (modes only), £ modcs , 
J2k n k UJ k (mean field), E m {. 



for the mean field alone and for the total two-point func- 
tion, using a log-log scale and leaving out the zero mode. 
While the high-momentum modes are still exponentially 
suppressed, the low-momentum modes have acquired a 
power law distribution. The quantum-modes-only dis- 
tribution does not behave as a power law (cf. Fig. |3J). 
The particle numbers as obtained from the mean field 
only and those including the modes have different pow- 
ers, — 1.5 and —0.67 respectively. Already much earlier, 
around tm = 8000, this distribution starts to emerge, 
with 20% larger powers. 

The power-law behaviour in the low momentum modes 
of the mean field apparently influences the quantum 
modes, in that their low momentum modes are enhanced 
in comparison to the classical T/tOk- We have seen this 
clearly in a plot of nkLUk {— * T for classical thermal 
equilibrium) showing a peak at k = and a "classi- 
cal plateau" at the interval k/m = 1.0... 2. 2. Simi- 
lar behaviour has also been found in the other runs at 
A/to 2 = 1/6, E/Lm 2 = 2 and A/to 2 = 1/8, E/Lm 2 = 4. 

In a purely classical simulation using the same set of 
parameters power-law behaviour is not found. This sug- 
gests that the interaction of the mean field with with 
the quantum modes plays a crucial role, even though the 
latter do not show power-law behaviour. 



2. Time scales for energy exchange 

In the previous section we obtained the scaling be- 
haviour cx A~ 3 for the time scale of approximate equi- 
libration based on the particle distribution. In this sec- 
tion we will use the energy density in the different parts of 
the field to find the short-time equilibration behaviour. 
Figs. 15171 show the results from one of the simulations 
at A/m 2 = 1/6, E/Lm 2 = 4, plotted at different time 
scales. In Fig. \5\ showing the early stage, the energy as 



+ N=128, Lm=32, X/m 2 =1/6, E/Lm 2 =4, 10 initial cond. 
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FIG. 6: Energy contributions for the same run as in Fig. 
for intermediate times. 
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FIG. 7: As in Fig.E|for long times. 
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obtained from the quasi-particles is also included. The 
quasi-particle picture appears to give a reasonable rep- 
resentation of the energies, with a roughly constant 10% 
mismatch in the total energy, to which we will come back, 
below. We furthermore see that the total energy in the 
quasi-particle picture is almost constant, corresponding 
to a quasi-particle number that is itself almost constant. 
This is consistent with the chemical potential found in 
the BE fits. We have checked that the dispersion rela- 
tion of the quasi-particle energies is close to that of free 
particles, u>k = \J m 2 s ± fc 2 , but with an effective mass 
m e ff that is larger than m, as can be seen, for example, 
from the minimum values of tOk m Fig.^ consistent with 
what follows from the effective potential. Coming back 
to the 10% mismatch in the quasi-particle energy: we 
have checked that it is neither a finite volume nor a finite 
lattice distance (spatial or temporal) effect. However, it 
does depend on the type of initial conditions used. For 
example a simulation with Bose-Einstein type initial con- 
ditions, as used in Ref. Q, with the same total energy, 
leads to a much smaller mismatch of almost 3%. It is 
important to note that we are looking at an effective de- 
scription of an interacting theory and it is not expected 
that the total energy in these particles is completely equal 
to the actual total energy as derived from the effective 
Hamiltonian. One therefore expects them to be closer 
when interactions are less important. As we will indeed 
see in Section UlI B 21 at smaller coupling and energy the 
mismatch becomes much smaller. 

Looking at the contributions E m f and -E mo dos to H c h, 
we see a relatively rapid transfer of energy from the mean 
field to the modes until a time of the order tm w 50. This 
exchange takes place fairly locally in momentum space, 
as is found by examining the mean field and mode contri- 
butions to rik, a phenomenon that we call local fc-space 
equilibration. At time tm ~ 100 most of the particle 
number already comes from the modes, whereas the to- 
tal distribution is still reasonably close to its initial form. 
After tm « 50, energy is still going to the modes, but 
with a slower rate. The behaviour in this second region, 
from tm 50 until tm 2000, (see Fig. |SJ) can be fitted 
reasonably well with an exponential form 

A + Be-^ T , (51) 

yielding Tin w 100 — 150. If we look at the long time 
behaviour, as plotted in Fig. [7| we see there is also a 
much longer time scale of the order 6000, on which energy 
is going back into the mean field. This time scale is 
comparable to the time scale of the emerging power-law 
behaviour, discussed in Section fill A II The appearance 
of this power law is accompanied by a large increase in the 
particle number in the zero mode of the mean field and 
therefore also in the average energy density of the mean 
field. We recall that classical behaviour only becomes 
visible at larger time scales of the order 15000. 

In order to make a quantitative comparison between 
different couplings and energies for the initial rapid ex- 
change of energy between mean field and modes, related 



rm E/Lm 2 = 1 E/Lm 2 = 2 E/Lm 2 = 4 
X/m 2 = 1/6 137 70 39 

{tm < 500) 

A/m 2 = 1/8 215 108 50 

{tm < 800) 

A/m 2 = 1/12 688 207 112 
{tm < 2500) | | | 

TABLE II: Initial energy-exchange time scales for the flat- 
ensemble initial conditions. 



rm 


Peak 1 & 2 


Peak 1 & 3 


"symmetric", Hartree 


160 ± 31 


360 ± 35 


"symmetric" , classical 


90 ±18 


156 ± 34 


"broken", Hartree 


49 ± 11 


84 ± 14 


"broken", classical 


41 ± 14 


63 ±16 



TABLE III: Auto-correlation times for flat ensemble type 
initial conditions. In all cases the coupling A/|/i 2 cn = 1/6. In 
the "symmetric phase" v 2 = 0, A/m 2 = 1/6 and E/Lm 2 = 4, 
whereas in the "broken phase" v 2 = 6, A/m 2 = 1/12, and 
E/Lm 2 = 0.5. 



to the local thermalization, we fitted the energy density 
in the mean field to a function of the form 1)510. The 
results are summarized in Table ITT1 Using the energy in 
the quasi-particle picture, instead of the effective Hamil- 
tonian, gives the same results. 

Leaving out the run at the lowest coupling and energy, 
X/m 2 — 1/12, E/Lm 2 = 1, which we did not see ther- 
malize, the time scale is roughly proportional to E~ l at 
constant coupling and to A -3 / 2 at constant energy den- 
sity: 

t- 1 ^Cm{E/Lm 2 ){X/m 2 ) 3/2 . (52) 

We have checked this behaviour explicitly by plotting the 
different energies as a function of [X/m 2 f/ 2 {E/Lm 2 )t 1 
which led to C = 0.10. 

For the lower energy density E/Lm 2 = 0.4 the results 
are very similar to our simulation of the Gaussian wave 
packet, which we describe in Section IlII Bl In particular, 
we encountered the local fc-space equilibration. If we ini- 
tially excite only a few modes this process can be seen 
even more clearly. We have not simulated long enough at 
this low energy density using the flat initial distribution 
to see the emergence of classical behaviour. 

3. Auto-correlation time scales 

To further investigate time scales, we also analysed 
the time-dependent auto-correlation function of the mean 
field, as in jlCJ. lllj . Using flat ensemble initial conditions 
in both the "symmetric" and the "broken phase" , with 
either Hartree or classical dynamics the auto-correlation 
function was obtained from the average mean field only: 
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N=128, l_m=32, X/m 2 =1/6, E/Lm 2 =4, Hartree, t =35000.. 70000 




FIG. 8: Auto-correlation function for the average mean field 
using the flat Hartree ensemble. The lower "arcs" are an effect 
of sampling the oscillations cx log | cos m e gt | at discrete times. 



is a much more striking difference between the "bro- 
ken" and "symmetric phase" results. At an 8 times 
larger energy, the auto-correlation time in the "symmet- 
ric phase" is not smaller, but instead larger by a fac- 
tor of 3 — 4. One would expect qualitatively the op- 
posite effect. For example, for a thermalized system at 
a temperature T the damping rate may be expected to 
scale, in the classical approximation, as (AT) 1 / 3 , and 
bluntly using the values AT/m 3 = (1/6)(1/0.41) ("sym- 
metric", Table0 and (1/12) 1.1 ("broken", g) would 
give T" symm " /T«b r okcn" = 0.62 instead of the factor 3 — 4. 

Comparing with the time scale for energy exchange, we 
see that at high energy density the damping time is 4 — 9 
times larger than the energy-exchange time (cf. Table ITTT1 
and the upper-right entry in Table ITT|) . The systemat- 
ics of this is unclear to us: at the lowest energy density 
(and smaller coupling) we find on the contrary that the 
damping time is about half the time scale for energy ex- 
change (see also Sec. lIII B~3l for the Gaussian wave packet: 

Tdamp ~ 3500, T oxch « 7000). 



C(t) = (</>(to - t/2)4>(t Q + t/2))° - d.c. (53) 

Here <j>{t) denotes the spatially averaged mean field, the 
long overline includes averaging over a large time interval 
(which greatly reduces fluctuations), and d.c. stands for 
the disconnected piece. Fig. [8] shows an example in the 
large-time region, where the particle-number distribution 
has the behaviour shown in Fig. The time average was 
taken over the region tm — 35000 . . . 70000, and ten ini- 
tial configurations from the flat ensemble. We recall that 
the evident damping is seen also upon using only a sin- 
gle configuration , it is not caused by the average over 
initial conditions. The dip-like structure can be under- 
stood as being caused by interfering "twin peaks" in the 
spectral function [T^j. The damping time is quantified 
using a "fit" of the form exp(— t/r) through the first and 
second peaks. Using the third peak would give a roughly 
twice as large r. The results are given in TableHn where 
the errors are obtained with the jackknife method [l2T |. 
In the table a comparison is made with results from the 
classical approximation, and with results obtained in the 
"broken phase" . 

In the "broken phase" there is hardly any difference be- 
tween the classical and the Hartree results, although the 
Hartree result seems to indicate a slightly larger value. 
We recall that the particle distribution in this case ap- 
proximates the Bose-Einstein form reasonably well, and 
furthermore, that the damping time is within a factor of 
two of the analytically computed value using perturba- 
tive quantum field theory in the two-loop approximation 



In the "symmetric phase" the Hartree result is roughly 
twice the classical value. It is hard to interpret this in 
any detail as the distribution function in the Hartree 
case is so "unconventional" (cf. Fig. 0} and also the 
classical case is far from thermalized. However, there 



B. Gaussian wave packet 

In this section we focus on the initial condition speci- 
fied by the Gaussian wave packet © with A/m 2 =0.1, 
Am 2 = 2 and $ = 2.60106 (this value appeared in the 
preprint version of 0), which gives an energy E/m 2 = 
12.6. We used a volume Lm = 32, giving an energy den- 
sity E/Lm 2 = 0.394 which is practically equal to the 
smallest energy density 0.4 studied in the previous sec- 
tion with the flat ensemble. It is however still an order of 
magnitude larger than the highest energy densities stud- 
ied in Our lattice size in this case was N = 256 lattice 
points, and the temporal lattice distance ao = a/10, as 
before. We checked for finite volume and discretization 
effects by using different parameters and found that they 
do not influence the results discussed. 



1. Particle distribution function 

The initial Gaussian wave packet spreads and oscillates 
in the course of time and after t > L/2 the packet meets 
itself through the periodic volume. This can be seen from 
the plots of the mean field (f>(x), see Fig. 1 in which 
we have verified. 

The initial wave packet © represents a pure state, 
which can still be analyzed in terms of particle numbers 
and frequencies obtained from the two-point functions, as 
in Eqs. (|14|l . It is interesting to compare the so-obtained 
rik with the coarse-grained particle distribution at later 
times. If we assume free-field evolution we can calcu- 
late rik analytically, and it turns out that its average 
over (half) an oscillation period is time independent and 
close to the initial distribution, for large volumes. As 
derived in Appendix lAl this free-particle distribution for 
the Gaussian wave packet as initial condition is given by 
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FIG. 9: Time development of particle number. The order 
in the key corresponds to the curves from top to bottom at 
ui k /m = 2.5. 



FIG. 10: The power spectrum S k (t) - S k (0) and n k (t) - 
l/2w fc (t) = S k (t) - l/2w fc for tm = 180 ■ ■ ■ 200. The line 
represents a power oc ui~ :i touching the negative values of 
S k (t)-S k (0). 



ftr - ^ 2 ^ TTV2 ^ A .. ( 54) 

In Fig. ED we plot this free form together with the par- 
ticle numbers obtained in a simulation. 

We find it quite remarkable that already the ear- 
liest (time-averaged) distribution deviates significantly 
from the initial form Q54J1. A closer look shows that 
this deviation originates entirely from the first period 
(tm = . . . 2tt) . After that short time the distribution 
is almost stationary. Only after a time tm — O(10 5 ) do 
we see deviations arise. However, in the mean time there 
is an extensive exchange of energy between the modes 
and the mean field: initially all particle number and en- 
ergy is contained in the mean field, while in the later 
stage it is just the opposite. 

After tm sa 30000— 40000 classical behaviour, i.e. n k — > 
T/u>k, starts to emerge: the lower-momentum modes be- 
come under-occupied, while the higher modes become 
over-occupied. At no stage does the distribution resem- 
ble the Bose-Einstein form i|50|l . We recall that also with 
the fiat ensemble we did not see quantum thermalization 
at similarly low energy densities. 

Bettencourt et al. [5| studied the power spectrum of 
the subtracted two-point function S k (t) — S k (0) at times 
tm < 200. This appeared to show power behaviour 
~ k~ 3 to fc~ 4 , which was interpreted as evidence for 
the absence of BE-likc thermalization. As mentioned 
above we also see no BE thermalization at this low en- 
ergy density, but we find that the power behaviour is 
not without ambiguities. The aim of the subtraction in 
S k (t) — Sk(0) was to eliminate the vacuum contribution 
l/(2\/m 2 + k 2 ) from Sk- At large k this is a rather deli- 
cate procedure. For instance, a quasi-particle behaviour 
S k (t) = (n k (t) + 1/2)/ y/m(t)'- l + k 2 with a thermal-like 
mass m(t) that is expected to be larger than m in the 



"symmetric phase" , would give a negative result at large 
k, S k (t)~S k (0) sa ~[m{t) 2 ~m 2 ]/Ak 3 , where we neglected 
an assumed exponentially small n k (t). 

We would like to stress here the good features of the 
observables n k and uj k defined in Eqs. i|14|) . In Fie. 1 101 we 
have plotted S k (t) — l/2uj k {i) = nk(t)/2uk{t), as well as 
S k (t) — Sk(0), for the same parameters ($ = 1.838526, 
Am 2 = 2, Lm = 128, N = 1024) and in the same time 
regime as used in (We averaged over tm — 180 — 
200, which hardly affects n k {t) /2u>k{t) as it is practically 
constant.) The plot of Skit) ~ Sk(0) looks very similar to 
the ones shown in 5]- There is a lot of scatter at large u>k 
(ss k) and a more detailed analysis shows negative values 
interspersed with positive values (indicated separately for 
the power spectrum). On the other hand nk/2tOk shows 
less scatter and is mostly positive (only for uj k > 4 do 
negative values occur). However, note that the larger u) 
region could be affected by lattice artefacts. 



2. Energy densities and time scales 

To get an estimate of the time scales involved we com- 
pare the energy densities in the mean field, in the modes 
and in the total field for the Gaussian wave packet initial 
condition, as we did in Sec. IIII A~2l for the flat ensemble 
at higher energy densities. For short times these are plot- 
ted in Fig. 1111 together with the energy as derived from 
the quasi-particle picture l|17fl . For long times they are 
plotted in Fig. [H 

We see that the quasi-particle representation of the 
energies is in this case extremely good, there is hardly 
any visible difference with the exact energies based on 

Furthermore, for early times fFig. lll|) there is an oscil- 
latory behaviour with a period tm sa 130. Note that, due 
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FIG. 11: The different energies for the Gaussian wave packet 
initial conditions at short times. From top to bottom: energy 
from modes and mean field, energy from mean field, energy 
from modes. 



FIG. 13: The auto-correlation function as determined from 
the average mean field only. The result when using classical 
dynamics is shown in the insert on a linear scale, with an 
exponential fit to the Hartree result. 
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FIG. 12: The different energies for the Gaussian wave packet 
initial conditions at long times. 

to the periodic boundary conditions on the system with 
size Lm = 32, the Gaussian packet already "meets itself" 
after a time tm = 16, much shorter than this resonance 
time. The resonance is caused by the difference between 
the effective mass terms of the modes and mean field, 
w 2\<fr 2 . This mass difference has a small value, fluc- 
tuating around 0.030 — 0.050, corresponding to a period 
210 — 126, approximately the observed period. 

The rate at which energy flows to the modes can be 
seen at long times fFig. H2|) . The energy in the mean field 
in the interval tm < 60000 can be fitted reasonably well 
to an exponential function of the form 151(1 . yielding an 
equilibration time scale rm ~ 7000, roughly two orders 
of magnitude larger than what was found in the "broken 
phase" at similar energy densities and couplings. 

Using a sum of waves as initial condition at similar 
energy density shows the same kind of resonance in the 
energy exchange between mean field and modes, with 



period ss 170 and mass difference fluctuating around 
0.02 — 0.04. Averaging over the flat ensemble the os- 
cillations die out after tm ss 4000 — 5000. 



3. Time scales from the auto-correlation function 

We also evaluated the auto-correlation function for the 
Gaussian wave packet. Since we do not average over ini- 
tial conditions we cannot calculate a statistical error. We 
therefore averaged over two different time intervals, giv- 
ing some idea of the size of the statistical uncertainty. 
The result is plotted in Fig. ED At this low energy, the 
damping time is roughly half the energy equilibration 
time. The inset shows - on a linear scale - the result for 
classical dynamics, using identical initial conditions. The 
exponential curve is a fit to the Hartree result: classically 
there is no visible damping. 

It would be interesting to compare the result with the 
fiat initial ensemble at the same energy density. However, 
for these low energies we need to simulate for very long 
times, which is quite a numerical effort. We therefore 
only calculated the auto-correlation times for the faster 
evolving high energy runs discussed in Section IlII Al 



IV. SCATTERING 

In this section we will discuss scattering features of the 
Hartree approximation. First we take a fresh look at the 
possibility of scattering via inhomogeneous mean fields 
by considering an initial state of two localized particle 
wave packets in position space. Next we give a perturba- 
tive explanation of "local /c-space equilibration" , our nu- 
merical result that the modes appear to equilibrate with 
the mean field primarily when they have the same wave 
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number. This occurs especially at low energy density and 
weak coupling and it has the effect that the initial particle 
distribution in the mean field is taken over approximately 
by the modes. We end this section with a brief discussion 
of higher-order scattering and thermalization. 



A. Scattering of two wave packets 

For homogeneous mean fields the Hartree approxima- 
tion cannot describe true scattering in which the mo- 
menta of the particles in the initial and final states dif- 
fer. However, in the inhomogeneous case the particles 
can truly scatter off the mean field. It is interesting to 
take an intuitive look at this in position space. Consider 
an initial two-particle state described by wave packets 
"01,2: 

|Vi^> = &%]6%]|0>, & + M (55) 

k 

This is a non-Gaussian and pure state, for which 

C rcn (x,i;x,i) = |^ 1 (x,t)| 2 + |^ 2 (x,t)| 2 (56) 

where 

" (57) 



V(x,t)= »£/ k (x,t). 



k 

Linearizing the Hartree equations in the "broken phase" , 
writing tf> = v + (f)' and keeping terms linear in <f>' while 
treating |?/?| 2 as being of the same order as </>', gives 

(d t 2 -A + mV = -3A(|^i| 2 + |^| 2 ), 
(df - A + to 2 + 6Aw0')^i,2 = 0. (58) 

If the wave packets approach each other within a distance 
of order 1/m they will scatter. 

So the interaction of the quantum modes with the clas- 
sical modes of the inhomogeneous mean field does lead 
to indirect scattering. Note that this happens especially 
in the "broken phase" : in the "symmetric phase" v = 
and the backreaction of the mean field disturbance 4>' to 
the particle waves i^n is suppressed. 



with w 2 = m? + fc 2 . will then play the role of an 
external field, as it is not altered by the interaction. 

We will show in the following that for not too large cou- 
pling and energy the equation of motion for g a reduces 
to that of a driven harmonic oscillator. Making use of 
the corresponding scattering diagrams we then conclude 
that, approximately, the only momentum modes of g a 
that are excited are those also present in the mean field. 
Since we will focus on the initial behaviour, when the 
system is still far from equilibrium and there is not yet a 
temperature, we will only use zero-temperature pertur- 
bation theory. 

We can write out the effective Hamiltonian in terms of 
the classical fields cf>, g a and the "external field" /°. In 
the "symmetric phase" and to second order in <? Q , this 
leads to the following interaction terms and correspond- 
ing vertex factors 



1 



-X<p 



3A^ 2 5>a° 5 ;) 

a 
a 

6A5>(/°<£)3?(/^) 

a, 



6A (60a) 

3A (60b) 

3A (60c) 

3A (60d) 



whereas in the "broken phase", writing cj> = v + (j)', we 
also have the three-point interactions 



Xv(/)' 3 

ex. 

3A^'^| 5q | 2 



6Xv (61a) 
3Aw (61b) 

3Xv (61c) 



In a first approximation we neglect the back reaction 
on the mean field and assume it is just oscillating around 
its minimum as a superposition of waves: 



4>(x,o) 



Ai sin(ajKit) cos(KiX — ipi) 



(62) 



B. Local fc-space equilibration 

To give an analytic interpretation of the "local fc-space 
equilibration" it is useful to focus on various interaction 
terms in the effective Hamiltonian (11611 . Although de- 
rived from a quantum system, this Hamiltonian can also 
be seen as describing interacting classical fields <f> and f a . 
It will be convenient to split the modes f a in a free part 
and a perturbation: 

f a (x,t) = f°(x,t) + g a (x,t), (59a) 

„ik a x — iuj a t 



fa{x,t) = 



y/2lV a L 



(59b) 



where ipi are random phases and = */ m 2 + Kf . 

The exact Hartree dynamical equation for the mode 
perturbation g a (x) in terms of its Fourier transform g a ^ 
is given by 



(df + cu 2 k )g ak = -3A / dx (<P(x) 2 + C Ien (x)) 




(63) 



Neglecting for the moment the higher order terms con- 
taining C rcn and g a k in the integral, the x integration 
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FIG. 14: Tree level scattering diagrams involving a single 
perturbation mode g a . Solid lines denote (j>, a dotted line f° , 
and the dashed line denotes g a . Time runs from left to right. 



can be performed, resulting in a sum over plane waves. 
The equation is that of a driven harmonic oscillator 



■IQa t 



(64) 



which leads to resonances that grow linearly in time for 
By inserting the explicit form i|62|l into Eq. I|63|) 



LOT 



n 2 . 



we find for each pair Ki 
lations: 



K,, four different resonance re- 



(65) 



(with uncorrelated ±), while the x integration gives four 
different momentum relations: 



k = k a + i]iK l + r]iKj, 



(66) 



where 771,2 = ±1. 

These two relations describe energy-momentum con- 
servation in scattering processes involving a single 4-point 
vertex, the interaction l(60b|l . Only 2 — > 2 processes in- 
volving this vertex can conserve energy and momentum. 
Furthermore, in 1 + 1 dimensions (since all particles have 
the same mass) it follows that the pair of incoming mo- 
menta must be equal to the pair of outgoing momenta. 
From the energy relation 165(1 it then follows that there 
are three possible diagrams, drawn in Fig. 1141 creating 
a g a particle with momentum k. The momentum rela- 
tion l|66|l now gives us three possibilities, 



k = f] 2 Kj 
k = r\xKi 



k a = —T]iKi (67a) 
k a = -mKj (67b) 

Vl K i = -V2 K 3 ( 67c ) 



For the last possibility k = k a , the 4>(x) 2 term con- 
tributes a constant. These terms only give rise to a time- 
dependent mass shift between the modes and mean field. 
We can conclude that to leading order the only excited 
modes are given by 



9±K i ,±K j 



(68) 



Note that only when just one mean field mode is excited 
(i.e. when i max = 1) the modes will remain diagonal. 
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FIG. 15: Leading 2 — > 4 scattering diagrams creating a g% 
particle in the "symmetric phase" . The intermediate line rep- 
resents the retarded Green function. 



We will now look at the neglected terms. The renor- 
malized mode sum C Ten (x) is equal to 

C ien (x) = J2(f**(x)9«(x) + fS.(x)g*(x) + \g a (x)\- 



(69) 

As we just showed, in lowest order, g a is only non-zero 
for k a € {Ki} and therefore the only non-zero Fourier 
components of C len (x) are the same as those in <fi(x) 2 : 
k = ±Ki± Kj. Therefore including the first order result 
for C Ten (x) in Eq. (|63|l will not change the set of excited 
modes. Finally, including the last term in Eq. I|63|l using 
the first order result (jrj5|) we can also find its contribution. 
The x integration gives a 5k> ,k±Ki±K • The frequencies of 
the correction to g a are therefore of the form WfcLfQiKj 
and we find exactly the same relation as following from 
Eqs. and TO. 

The above treatment can be extended by making a 
systematic expansion in A, <f> = 4>q + \<f>\ + \ 2 4>2 + • • • , 
fa = f a + A/j, + \ 2 fa + • ■ • , and using Green function 
techniques along the lines of |l3[ . 

As a check we performed a simulation, exciting only 
two modes K\ and Ki at low energy (E/Lm 2 = 0.04) and 
small coupling (A/m 2 = 1/12). The assumption of a free 
oscillating mean field turned out to be extremely good. 
We also checked the explicit form of one of the modes by 
examining \ftci \ 2 ■ We expect / to contain the two Fourier 
modes K\ and K2, and therefore \f\ 2 to contain momenta 
2Ki, 2K2, K\ + K2 and K% — K2, which were indeed 
the only modes found. In similar simulations at higher 
energy we found that the back reaction to <f> became more 
important, but the set of excited modes remained the 
same. 



C. Higher order scattering 

In order for the system to thermalize it is necessary 
that particles can change their momenta by scattering. 
As mentioned above, 2 — > 2 scattering cannot change the 
initial momenta in 1+1 dimensions. (With re-summed off 
shell propagators this restriction does not apply, cf. the 
nontrivial spectral functions in }10| and |14j , and the ther- 
malization found in Q.) Therefore at least one extra ver- 
tex is needed. In the "symmetric phase" only four-point 
vertices exist, as in Eqs. (I6U|I . The interaction i|60b|) is 
leading over l|6Uc|l and (|60d|) . because it is first order in 
g a , and since initially all energy is in the mean field, the 
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FIG. 16: Leading 2^3 scattering diagrams creating a g% 
particle in the "broken phase" . 

leading contribution to g-particle production comes from 
the two diagrams in Fig. El 

At this point it is interesting to realize what happens 
if the mean field is homogeneous. In that case g a always 
carries momentum k a . For inhomogeneous systems this 
restriction is lifted and thermalization becomes possible. 

In the "broken phase" both the couplings l|60[l and i|61|) 
contribute and there are three- and four-point interac- 
tions. The leading contribution to (/-particle production 
in this case comes from the two diagrams in Fig. ^j] In- 
tuitively one expects the finite range of the interaction in 
the "broken phase", due to off-shell particle exchange, 
to lead to more efficient thermalization than the zero 
range interaction in the "symmetric phase" . This is in- 
deed what we observed. 



V. DISCUSSION 

We will start this discussion with a summary of the 
behaviour at high and low energy density. This appears 
to be the distinguishing criterion for the thermalization 
behaviour of the Hartree approximation for inhomoge- 
neous systems, rather than, for example, an initial state 
being pure, as for the Gaussian wave packet, or mixed, 
as for the flat ensemble. 

At high energy density E/Lm 2 > 1 we see that the 
distribution rik acquires features of a thermal quantum, 
i.e. Bose-Einstein distribution. There is a time- and 
coupling-dependent chemical potential, of order unity in 
mass units. The temperature is roughly proportional to 
\J E I L and independent of the coupling. The coupling 
determines the time scale on which the approximate ther- 
malization becomes visible. The initial rapid exchange of 
energy between modes and mean field occurs on a time 
scale described by Eq. I|52|l . The quasi-particle picture is 
reasonable and the total particle number is very constant, 
in correspondence with the chemical potential. However, 
there is a mismatch between the total energy as derived 
from the effective Hamiltonian and that obtained from 
the quasi-particles, Eqs. (|16fl and (|17|) . As we already 
pointed out, this mismatch depends on the initial condi- 
tions, for example a Bose-Einstein type initial conditions 
cf. Ref. leads to a much smaller mismatch. For an 
interacting theory it is not expected that a quasi parti- 
cle picture based on free particles would give precisely 
the same value as what follows from effective Hamilto- 
nian and indeed the mismatch was practically invisible at 



lower energy and coupling. Furthermore, a considerable 
amount of the total energy comes from the zero mode, 
about 10%, making it very sensitive to the precise par- 
ticle distribution and probably making deviations from 
the quasi-particle picture more pronounced. 

After a very long time the energy flows back into the 
mean field, accompanied by the emergence of a power-law 
distribution for rik as a function of momentum k. Such 
power-law behaviour has not been found in the "broken 
phase" , and also classical simulations do not show such a 
behaviour, indicating that the back reaction of the modes 
plays an essential role. It is interesting to note that Boy- 
anovsky et al. ^| also found power-law behaviour for the 
occupation numbers, although it is unclear if the same 
mechanism is behind their finding. In their study of a 
3 + 1 dimensional O(N) model in the large N approx- 
imation, the power-law is caused by a nonlinear reso- 
nance of the back reaction of the modes on themselves, 
with terms of the form l/(u>k — m), diverging as 1/fc 2 in 
the limit k — > 0. This would give a 1/fc 4 behaviour for 
the particle number, different from what is found here. 
Furthermore, we only find power-law behaviour in the 
total- and mean-field-particle numbers, but not in that 
of the modes. The difference in power and the absence 
of power-law behaviour in the modes makes it improba- 
ble that the physical mechanism behind the resonances 
is the same. 

At low energy density E/Lm 2 <C 1, for the flat ensem- 
ble as well as for the pure-state wave packet, we do not 
find approximate thermalization to a BE distribution. 
Instead, the form of the total distribution rik remains 
the same for times that are many tens of thousands in 
units of m . The distribution then slowly turns over 
into a classical distribution. However, there is still an 
extensive exchange of energy between the mean field and 
the modes, leading to what we call local fc-space equili- 
bration, for which we were able to find an interpretation 
based on a perturbative calculation. At low energy, the 
time scale for this process is much longer than would fol- 
low from Eq. l|52|l . found at high energy densities. Fur- 
thermore, at short times the energy densities in the mean 
field and modes separately show a remarkable oscillatory 
behaviour, not seen at higher energies, which is caused 
by a difference in the effective mass of the mean field and 
modes. 

We obtained several time scales: for approximate Bose- 
Einstein thermalization, for the early-time exchange of 
energy between mean field and modes, for the auto- 
correlation function and for the evolution to a classical 
distribution. Most of these are much longer in the "sym- 
metric phase" than in the "broken phase", due to the 
absence of the finite range interaction. 

According to the BE-thermalization-time scale in 
the "broken phase" is of the order 25 — 35 for E/Lm 2 = 
0.5, while here, in the "symmetric phase", totbe = 
1500 - 1600 for E/Lm 2 = 4 (both at A/|/4J = 1/6). 

The energy-exchange time scale in the "broken phase" 
gives a result that is close to tbe, whereas in the "sym- 
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metric phase" it is related to local A:-space equilibration, 
much shorter than tbe> and it shows the behaviour (|52H . 
For E/Lm 2 = 0.5 and A/|/4J = 1/6, Eq. gives 
Tin ps 300, much longer than the 25 — 40 we found in the 
"broken phase" at the same energy and coupling. 

Also the damping time, obtained from auto-correlation 
functions, is much longer than in the "broken phase" , 
even at much higher energy and larger coupling. Com- 
pared to the value obtained using classical dynamics, it 
is roughly twice as large. In the "broken phase", both 
values are comparable in size. At low energies the damp- 
ing time seems to be much longer, but this needs more 
study. 

The last time scale is that of classical equilibration. 
Since we are just solving a large number (2N 2 + 1) of 
local classical non-linear equations, one may expect clas- 
sical equipartition to set in at some point. This equiparti- 
tion is, however, non-trivial because of the large number 
of conserved charges Q • For example, the emerging clas- 
sical temperature is of order E/N and not E/N 2 Q- 3 
Depending on energy and coupling we can already see 
a first emergence of classicality at times rm = O(10 4 ). 
This is still about an order of magnitude longer than 
what was found in the "broken phase" in |3j. However, 
full classical equilibrium is expected only for huge times, 
much larger than the rm — O(10 6 ) found in the "broken 
phase" for an artificially small system at E/Lm 2 = 36 |2j, 
and beyond the already large times of order 10 5 reached 
in this study. 

Remarkably, the equilibration time scale found in [l| 
using classical dynamics appears to be shorter. The em- 
pirical formula [l| 

l/mr class = 5.8 10- 6 (6AT/to 3 ) 139 , (70) 

with T = E/N the classical equilibrium temperature, 
would give equilibration times tm = O(10 5 ) - O(10 7 ) for 
the various parameters used here. This difference in time 
scales can be interpreted as follows. Classical dynamics 
has also be studied in the Hartree approximation, and 
the latter shows up as an unstable fixed point of the 
full dynamics [lj. This Hartree fixed point depends on 
the initial conditions. In our case the mode functions 
are initialized with quantum- vacuum form and the 
resulting dynamics (seen as a classical system with order 
N 2 fields) appears to linger for a very long time near 
a Hartree fixed point, longer than when using classical 
dynamics. 

For our inhomogeneous initial conditions we have not 
been able to pin down the fixed point analytically, but in- 
tuitively one may expect the system to be close to it when 
the mean field has lost most of its energy and has started 
fluctuating about a homogeneous average. Making a ho- 
mogeneous approximation to this situation would lead 



to a Hartree stationary state. Such a state can have an 
arbitrary particle distribution rik, which, given our out- 
of-equilibrium initial conditions, turns out to have BE 
features when the energy density E/Lm 2 3> 1. Appar- 
ently, when the energy density is small, E/Lm 2 <C 1, the 
system leaves the fixed-point region before BE-like ther- 
malization sets in, because we have seen only classical- 
like equilibration emerging in this case. 

Finally, we comment on the results of Bettencourt et 
al. |5J. As mentioned in Sec. IIII 5T1 we have essentially 
confirmed their numerical results. The energy density in 
the simulations in [I| was rather low, namely E/Lm 2 = 
0.00042 and 0.0045, so in view of our results summarized 
above, no sign of a BE distribution is to be expected 
with the Hartree approximation at the times tm < 200 
covered in |5j, nor at any time later. 

It is remarkable that we do find Bose-Einstein be- 
haviour at larger energies E/Lm 2 ^> 1, but of course, the 
fact remains that one needs to improve on the Hartree 
approximation in order to achieve thermalization at all 
energies. This may take huge times at low energy densi- 
ties. 

It has been remarked [f| that the Hartree approxima- 
tion is expected to be valid up to times tm ~ m 2 /\ = 
O(10). We agree with this statement when applied to 
the detailed time-evolution of observables, but it does 
not necessarily apply to observables such as our quasi- 
particle distribution n/-(t) or energy tOk(t), which are 
coarse grained in time and space and/or averaged over 
initial conditions in the Hartree ensemble approximation. 
For comparison, consider a gas of classical point particles 
with Lennard- Jones interactions. Any numerical approx- 
imation to the detailed time evolution will soon go dis- 
mally wrong due to the chaotic nature of the system, but 
this does not preclude an accurate evaluation of, say, a 
coarse-grained particle-distribution function. With this 
in mind we have studied our system for times as large 
as seemed necessary, which led to very large times in- 
deed. First experience 0] indicates that the situation 
is not very different in 3+1 dimensions, where also large 
equilibration times may be expected for the c/> 4 model at 
moderate couplings and energy densities. 
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APPENDIX A: PARTICLE NUMBER OF THE 
GAUSSIAN WAVE PACKET 



J For N spatial lattice sites there are 2N 2 real de grees of freedom 
in the mode functions. 



We calculate here the initial two-point functions for 
the Gaussian wave packet initial condition @, the cor- 
responding particle number rik and energy Uk , and their 
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subsequent free field expressions. The calculations will 
be made in the continuum limit, in our finite periodic 
volume. 

The mean field contributions to the two-point func- 
tions are given by 



S(x,y) mt = 4>(x)cj>(y) - <p(x) cf>(y), (Ala) 
£/(x,y) mf = 7r(x)7r(»)-^(5) (Alb) 



where at first we shall average only over space, i.e. 



1 



x)<p{y) = 7" / dz (f>(x + z)(p(y + z). 



(A2) 



The initial mean field is given by Eq. JjJJ, or in terms of 
its Fourier transform: 



dx e~ lkx (j){x) = $V2^I e~ k A l 2 . (A3) 



Adding the contributions in Eqs. (|A5|) and (|A6|I and ap- 
plying the definition (|14|l . the initial instantaneous par- 
ticle number and frequency become 



nfc(O) 
Wfc(O) 



1 - 1 



,(o) 



(A7a) 
(A7b) 



Using free field dynamics the instantaneous particle num- 
ber would get an oscillating component according to 
Eqs. i|A5(l . If we also course grain in time, the discon- 
nected parts of S and U vanish, while both cos 2 and 
sin 2 1/2. We then find 



Since ir(0) = 0, the free-field (i.e. for A — > and /i 
Mrcn = m) evolution of 4>k is given by: 



<t>k(t) = 0fe(O)cos(4 O) i): 



(A4) 



where to^ — \J m 2 + k 2 . A straightforward calculation 
gives 



sT f = (l - <5 fe ,o) 



2 cos 2 (4 o) t) 



7 fc 



(4 0) ) 2 ^sin 2 (4 ^) 



(A5a) 



(A5b) 



_,frec 



2L 



J k i 



(A8) 



which are time independent. 



For large volumes 2 u^ (j) 2 k /L <C 1, expressions (|A7a|) 
and i|A7b(> reduce to (|A8|) . For the parameters as used 
in SectionHTH A = 2, $ = 2.60106, Lm = 32, we have 



9 (0) ,2 

« 5.30 + k 2 /m 2 e - 2k2 ' m \ (A9) 



where the delta functions come from the disconnected 
pieces. The modes just contribute the vacuum fluctua- 
tions: 



^finodes 



1 



2uo 



(0)' 



jo) 

modes k 



(A6) 



Plotting rifc(O) (or log(l + l/nfe(0)) versus cjfe(0) we find 
that this only compares well with a similar plot of n^ 00 
versus for fc ^ 2m. So at times tm S> 1 it is best 
to use the time-averaged free-field determinations for the 
comparison with the interacting Hartree evolution. 
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